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We examine the question of whether fluids and crystals are differentiated on the basis of their 
zero frequency shear moduli or their limiting zero frequency shear viscosity. We show that while 
fluids, in contrast with crystals, do have a zero value for their shear modulus, in contradiction to 
a widespread presumption, a crystal does not have an infinite or exceedingly large value for its 
limiting zero frequency shear viscosity. In fact, while the limiting shear viscosity of a crystal is 
much larger than that of the liquid from which it is formed, its viscosity is much less than that of 
the corresponding glass that may form assuming the liquid is a good enough glass former. 

I. INTRODUCTION 

Response theory and Green-Kubo relations provide a good understanding of the microscopic origins of viscoelasticity 
in fluids^^. They show that all fluids are in fact viscoelastic. However the range of frequencies over which one sees 
a crossover from low frequency viscous behaviour to high frequency elastic behaviour, varies by more than 10 decades 
for various common fluids. Solids are also expected to be viscoelastic, exhibiting viscous as well as elastic behaviour—. 
Of course there is a limit to the maximum strain amplitude that can be applied to a crystal before it fractures, 
cleaves, plastically deforms or melts^. However provided this limit is obeyed one can in principle compute or measure 
elastic constants at all frequencies including zero and shear viscosities at all frequencies except zero. This presents an 
interesting question: how does the limiting zero frequency shear viscosity of a crystal compare to that of a fluid or a 
glass? 

In contrast to fluids the microscopic origins of the rheological properties of crystals are poorly understood. There 
is no equivalent to Green-Kubo theory that has been proposed for the viscoelastic behaviour of crystals. We know of 
no collection of the frequency dependent shear viscosity of crystalline materials. There are collections of experimental 
data for the frequency dependent shear moduli. 

It is easy to understand why this surprising state of affairs exists. The standard derivation of Green-Kubo ex- 
pressions for the Navier-Stokes transport coefficients relies on the Onsager regression hypothesis and a solution of 
the fluctuating Navier-Stokes hydrodynamic equations^. For crystals the corresponding elasto-hydrodynamic mode- 
coupling equations are quite complex (and anisotropic). As far as we are aware no one has derived the Green-Kubo 
relations for the frequency dependent elasto-hydrodynamic coefficients of a crystal. Furthermore when you shear a 
crystal the underlying equilibrium state varies with the strain. This is not so in a fluid. This difference implies that 
the Green-Kubo derivations themselves are inherently more complex for crystals than for the corresponding liquid or 
gas. 

We adopt a much easier approach to solve this problem. We restrict ourselves to the case where there is no significant 
linear creep and then employ the SLLOD equations for time dependent planar Couette flo\4. These equations give an 
exact description of adiabatic time dependent planar Couette flow arbitrarily far from equilibrium. These equations 
convert a technically complex thermal transport process into a much simpler mechanical process that can be analysed 
using a thermostatted version of Kubo's response theory^. 

II. THEORY 

A. Thermostatted SLLOD equations for isothermal planar shear 

We use the standard isokinetic equations of motion which feature a synthetic thermostat under the condition that 
the total peculiar momentum is always zero. We know some important facts about synthetic thermostats of this 
type: 

1) Their equilibrium distribution function is known&il. 

2) Any artifacts in the dynamical correlation functions due to the synthetic thermostat are at most of order 
0(1 /N) where N is the number of particles in the system^. 

3) The linear response of the system to an external field is devoid of artefacts due to the synthetic thermostat. 



4) If we wish to study nonequilibrium phenomena outside the linear response regime we can arrange things 

such that the thermostat only acts on a region far removed from the system of interest. 

The isokinetic equations of motion are 

q, = pi/mi + fi< (T) -F e 

p 4 = F 4 (q)+D, t (r)-F e ~a Pj 

EjIiPi-Fi + Pi-D i (r)-F e 

a = ^jv w 

Ei=i p* ■ Pi 

where and Pi are the position and peculiar momentum of the ith particle, q is the 3 AT dimensional vector of all 
the positions, T = (q, p) is the 6N dimensional phase space vector, m; is the mass of the i th particle, F, is the force 
on the i th particle due to interactions with other particles, F e is the external field which drives the system away 
from equilibrium. When it is set to zero, if the system is ergodic and has decaying memory, any arbitrary initial 
distribution will eventually relax to equilibrium^. Cj and Dj are second rank tensors which couple the system to the 
external field and a is the thermostat multiplier which holds the value of the peculiar kinetic energy, YliLiPi/^m, 
constant. The equilibrium distribution function (when F e = 0) for the system is given by 

/(r) = S(K(T) - Ko)S(p M ) exp[-/3ffo(r)] 

Z = J dT5(K(T) - K )5(p M )eM-(3Ho{r)] (2) 

where K(F) is the peculiar kinetic energy which is fixed to the value Kq, Hq(T) — K(p) + 3>(q) is the Hamiltonian 
with $(q) the potential energy due to the particle interactions, pm — Ei=i Pi ls the total peculiar momentum, 
and /3 _1 = fcsT = 2Kq/(3N — 4) where fcs is Boltzmann's constant and T is the equilibrium thermodynamic 
temperature. We note that if the system is ergodic with decaying memory the distribution function Eq. <[2j) is the 
unique, dissipationless equilibrium stated. We also note that Eq. |f2|) includes all finite size corrections. The Helmholtz 
free energy of the system is£ 

A=-k B T\n[Z]. (3) 

For the system undergoing planar shear we use the so called isokinetic SLLOD equations of motion^, 

q, = Pi/mi + ij(t)q yi 
Pi = Fi - ij(t)p yi - ap, 

n _ Sili Pi ■ F i - i(t) PxiPyi 
S»=i p» • p* 

where i is the unit vector in the direction of the x Cartesian axis, p X i is the x component of Pi the peculiar momentum 
of particle i, and j(t) — du x /dy is the time dependent strain rate where u x is the x-component of the steaming 
velocity. In computer simulations Eqs. Q are used in conjunction with Lees-Edwards shearing periodic boundaries^ 
to minimize the system size dependence of the results. The adiabatic form of these equations give an exact description 
of adiabatic shear flow arbitrarily far from equilibrium. They are equivalent to Newton's equations of motion plus 
an integrated shift in the x-laboratory velocity of J Q ds ; y(s)q y i(s) for every particle. Decomposing the strain rate into 
an infinite sum of infinitesimal Heaviside steps shows that the adiabatic form of Eq. (0]) is exact for time dependent 
planar Couette flow. 



B. The difference between solids and fluids under quasistatic strain 

A fundamental difference between a fluid and a solid is that while a solid can support a small externally applied 
stress indefinitely, a fluid cannot. Fluids will always flow in response to the applied stress thereby eventually reducing 
the magnitude of the stress to zero. If we subject a liquid to a quasistatic strain rate there will be no work done in 
shearing it. To prove this we note that the rate at which work is done in shearing a single ensemble member is given 

by 

W(t) = HS d (T) = ^.f» d ^ ~J(r)VF e (i) = -j(t)VP xy (T), (5) 



where the flux, J, introduced here is defined by Eqs. |(4]&[5|) and P xy is the xy element of the pressure tensor (closely 
related to the shear stress, a xy — — (P xy )). Using Eq. we see that 

N N 

r/ , v;| , ; -Y. I' (6) 

1=1 1 i=l 

To leading order we have (P xy ) ~ 7 + C(j 3 ) and so if we calculate the work required to quasistatically strain a fluid 
by a fixed amount J7, at constant shear rate 7, we obtain 

{AW) QS = - lim W / dt (P xy (t)) = Urn VtffSi = 0. (7) 

where 77 is the shear viscosity. This is relevant because the change in free energy due to the strain is exactly given by 
the quasistatic work done and thus the underlying equilibrium free energy of a fluid does not change with a strain. 

Let us now consider what happens if we subject an initially unstressed solid, 7 = 0, to an infinitesimal change in 
strain. We assume that the final strain £7 is sufficiently small that the solid responds according to linear elasticity 
theory. Thus the average shear stress is related to the zero frequency shear modulus, Go, and the strain by the 
equation (P xy (t)) = — Go"f(t). The change 5j may be effected by perturbing the boundary conditions, 



/•<5-y/7 

{AW) QS = - lim W / dt (P xy {t)) 
7^0 Jo 



r&lh 1 

= limVGo7 2 / tdt = -VG 6j 2 . (8) 
7^0 Jo 2 

Because a solid can support a stress for an indefinite time the underlying equilibrium free energy will depend upon the 
change in strain, as in turn will the partition function, Eq. |(3]). Thus the expression for the equilibrium distribution 
function will now explicitly depend upon the strain through the partition function, Z(Sj), and so will the equilibrium 
average of a phase variable, (^(r))^ . 

We wish to calculate (to leading order) the change in the xy-element of the pressure tensor for a solid subject to a 
shearing deformation with a strain £7. The equilibrium average can be calculated from the expression, 

, p . J D{67) dTP xy (T)ex P [-l3(H (T)} 

^*r«- J D(Si) dTeM-(3(H (T)j ■ (9j 

In this equation D(Sj) defines a phase space domain which is strained an amount £7, from a reference domain D(0). 
This may represent changes in the boundary conditions. The average (. . .) eq is an equilibrium average taken over 
the domain D(0) and the average (■■•)s~ eq is similar but taken over the domain D(Sj). Because of the anisotropy of 
crystals the observed stresses will be strong functions of the alignment of the crystal relative to the strain or strain 
rate, tensor. For simplicity we do not use notation that makes this alignment explicit. 
The transformation between the two domains is given by the equation, 



r' = r - sr (10) 



where ST is 



ST = S~/(q yl ,0, 0, q y2 , 0, 0, . . . , q yN , 0, 0, 0, 0, 0, . . . , 0, 0, 0). (11) 
We can transform the average Eq. ^ using the coordinate transformation as 



f dT' 

JD(0) 


ar 
ar' 


P 




' + 5T)exp[-/3(H (T' + 5T)] 


f dT' 

JD(0) 


ar 

ar' 


exp[-/3(fl (r / + «r)] 



(P \ - __w i__ 1 g ' : (12) 

Noting that the Jacobian is unity, that dT' is a dummy integration variable and expanding P xy and Hq to leading 
orders in 67 gives 

f D(0) dT[p xy (T) + ■ \7p xy (T)} ex P [-/3([g (r) + ■ yg (r)]] 

{ xv) *r><*- f D{0) dTexp[-l3([Ho(T) + Sr-VH Q (r)]] ' (13) 



Approximating the exponentials to leading order in ST gives, 

_ f m dT [P xy (T) + ST ■ VP xy (T)} [1 - (38T ■ VH (T)} exp[-/3(ff (r)] 



{Pxy) 



V I S-y,eq 



J D{Q} dTexp[-(3(H (T)} - f m dTST ■ VH (T) exp[-/3(iJ (r)]] ' 

and expansion of the denominator, to leading order in ST, gives, 

J D(Q) dT [P xy (T) + ST ■ VP xy (T)] [1 - (3ST ■ VH (T)} exp[-(3(H (T)] 



{Pxy) 



Sd(o) dTexp[-P(H (T) 



(14) 



x 1 + 



0f D(o) dTST • VH (T)ex P [-(3(H (T)] ' 
J D{0) dTexp[-P(H (T)] 



Applying the coordinate transformation, Eqs. lfT0|) & (fTTj). we see that 

ST • VH (T) = -S 7 P* y (T)V, 
where P xy is the configurational component of the xy-element of the pressure tensor and 

N N dF x 



ST ■ VP xy (T)V = ^EE 777" '/«•'/.•/ .•• 



i=l j=l 



Ixj 



Substitution into Eq. (fl5| gives 



(15) 



(16) 



(17) 



+^V 2 [(P xv P* y ) 



0,eq 



r / \ ~* \ 9F x i 



{Pxy) Q eq (Pxy) eq 

> + O (<5 7 2 ) 



> l=1 J = 1 ' 0,eq 

,1 l-n \2 



= (P*y)o,e q V + ^iV 2 [(P 2 xy ) eq - (Pxy)i eq \ - S 7 V ( 9oo )^ eq + O (Sj 2 
Where for simplicity we define the phase function by the equation 



N 2 
Pv 



N N 



^— ' rrii ^— ' ^— ' 

i— 1 i—1 j=l 



dF x 



-QyiQyj 



lxj 



which implies the zero frequency shear modulus is 



Go - (<U ,e„ ~ PV 



(P*v)o,, 



0,eq 



(18) 



(19) 



(20) 



The zero frequency shear modulus is thus the sum of a fluctuation and a nonfluctuating term. The nonfluctuating 
term (valid at zero temperature) was given by Born 9, lfi in 1939. The derivation of the correct finite temperature 
result was first given by Squire et. al . 11 ' 12 in 1969 and rediscovered in 1986, see ref. [H. For a fluid, the sum of these 
two terms is exactly zero since the shear modulus is zeroi^. For a solid these two terms do not cancel and there is a 
non-zero shear modulus. 

To gain a better understanding of these two terms consider the response of a system to an impulsive strain rate: 
j(t) — S~/S(t). From the SLLOD equations of motion Eq. (j4]) we see that for impulsive shear the change in the phase 
space vector is 



ST = 5j(q y i,0,0,q y 2,0,0, . . . , q V N,o,o, -p y i,0,0, ~P V 2,0,0, . . . , -p yN ,0,Q). 



(21) 



If we now consider some phase variable B(T) whose functional form is not explicitly dependent on the strain, 

B(r(0+)) = B(T(Q-)) + VB(T(Q-)) ■ ST(Sj) + 0(5j 2 ), (22) 



and substitute P xy V for B we see that 



(^(0+)) V = (P xy ) 0eq V - 6jV ( 9oo ) 0<eq = (P xy ) Qeq V - SjG^V, 



(23) 



where (P xy (t)) is the nonequilibrium average taken, at time t, in this case at time t = + which is directly after the 
system has been subjected to the impulse. The nonfluctuating component of the zero frequency shear modulus is in 



fact the infinite frequency shear modulus, Goo, so that, Gq = G c 



(3V 



. For all systems 



0,eg 



(solids or fluids) the infinite frequency shear modulus is given by Eq. lfT9]). At infinite frequency there is no time 
for the system to recognize whether it is a fluid or a solid. In a fluid, and only in a fluid, the zero frequency shear 
modulus is zero and hence we get a second exact expression for the infinite frequency shear modulus that is only valid 



for fluids: G^= (3V{ P xy - (P xy ) 



'0,eq 



Here the superscript F indicates that this expression is only valid 



0,eg 



for fluids. This latter expression is familiar to those acquainted with the Green-Kubo expressions for the frequency 
dependent shear viscosity of fluids. 

Thus regardless of what state of matter we are dealing with 



0V 



p 



(P 



x V'0,eq 



G c 



G Q . 



(24) 



0,eq 



However the difference between a solid and a fluid is whether the zero frequency modulus is zero as in a fluid or a 
positive number as in a solid. 

It is important to remember that this perfect mathematical cancellation between the fluctuation term and the 
nonfluctuating term in fluids is no simple mathematical identity. Consider two systems with identical Hamiltonians, 
densities and equations of motion. The only difference between the two systems is their temperature. One is in the 
solid state phase and the other is in the liquid phase. You can transform between the two states by simply changing 
the temperature. Yet in the liquid the cancellation is perfect whereas in the solid it is not. This cancellation is a 
symmetry that is particular to the fluid state. 

The standard derivations of linear response theory assume that the underlying equilibrium distribution function 
does not change with strain. Clearly if we wish to treat a solid phase we must account for the effect of the underlying 
equilibrium distribution function depending on the strain. 



C. Linear response to shear for systems that are initially at equilibrium 

In this section we will consider the linear change in the stress in response to a small applied strain, firstly for an 
ordinary fluid, then a crystal which is initially unstrained and finally for a crystal with an initial strain that is not 
zero. 

Linear response theory is described in terms of a field and a conjugate flux. For the special case where both the 
field and the flux, which are vectors, have a common direction, only their magnitudes are relevant and we may define 
the flux, J, as, 

where V is the system volume, F e is the magnitude of the field which appears in Eqs. HJ and Q is the rate at which 
heat is exchanged with the synthetic thermostat. To leading order in the field (q(T) — Hq(T)^ = 0{F%) and thus as 
the field approaches zero so does the flux, J(T) = 0(F e ). 



1. Equilibrium fluid 



Consider a fluid in the isokinetic ensemble, Eq. ([2]), which is initially in equilibrium and then perturbed by an 
external field at time t — 0. Linear response theory gives,—, 

(B(t)) = (B) eq - (3V ( ds (J(-s)B(0)) eq F e {t - s), (26) 
Jo 

where B is some arbitrary phase variable, ( J{— s)B(0)) = {J(T(~s))B(T)) eq and r(— s) is the point in phase space, 
such that if we start at it and run the equations of motion forward in time (with F e = because the average is an 
equilibrium one) we arrive at the point T at time 0. Here we are interested in planar shear with B(T) = J(T) = P xy (T), 
{Pxy) eq — and F e (t) — j(t). So we obtain 

(P xy (t)) = -PV f ds (P xy (- s )P xy (0)) j(t-s). (27) 
Jo 



2. Initially Unstrained crystal 

For the case of a crystal, the underlying free energy changes with the strain due to the change in the boundary 
conditions and a stress can be supported indefinitely. Of course there are other processes than just planar shear which 
could result in this behaviour and so we will represent the change in the free energy using the arbitrary parameter A. 
In the case of planar shear we will have A = 7. We have recently given a generalisation of linear response theory for 
such a case where a system may be simultaneously subject to a dissipative field, F e , and a parametric change, A(i), 
to its equilibrium stated. In this paper we proved that to linear order in an arbitrary dissipative field and parameter 
the average linear response of a phase variable B(T, A) that may depend on the parameter is, 



(B(t)) m = (B) m , eq -PV ds (J(-s)B(0)) HOUq F e (t-s) 



Jo 



(3 / ds 



dA(\(0)) / dH(T(-s),X(0) 

— {B(T, A(s))) A(0) eg - ( — B(T, A(s)) 



\(0),eq 



X(t-s). (28) 



We choose to set the parameter as the strain and the dissipative field as the strain rate. In this case the Hamiltonian, 
H , and the phase function B(T), have no explicit dependence on the strain or strain rate, but averages will still be 
dependent on these parameters via the boundary conditions. Applying Eq ([28|) to this situation gives 



<i?(t)) 7(t) = (B)^ t)>eq -0V j n ds{P xy (-s)B(0)) i{0) eq j(t-s) 

,&4(7(0)) 
(?7 



~P Jl ( B )j(0),e Q I dsi(t-s), (29) 



where (B(t))^^ is the value of the nonequilibrium average at time t. 

Let us now consider the example of a planar shear impulse again, and choose B(T) = J(T) = P xy (T), j(t) = V t < 
and j(t) = 7i<5(i). Because the initial stress is zero we will have (P X y) 7 / t \ eq = V t < 0. We consider the terms on 
the right hand side of Eq. lj29j) in turn. The first term is easily seen to be 

<^) 7l , eg = -Go7i + 0(7i 3 )- (30) 

The second term is easily evaluated as —0Vji (P xy (T(—t))P xy (T)) Lastly we need to calculate the change in the 
equilibrium free energy caused by the strain. We note that Eq. ([3]) relates the free energy to the partition function. 
The change in the partition function is easily computed as 

Z(N,V,T,~/ 1 )= [ dTexp[-p(H (T)]-p [ dT6T ■ VH Q {T) exp[-/3(F (r)]] 
Jd(o) Jd{o) 

= Z(N,V,T, 7 (Q-)) [ 1 + ^7i(^*> 7(0 - ) , e J , (31) 



and the free energy is thus 



A(N, V, T,7i) = A(N, V, T, 7 ((T)) - lx V <P* ) Q . (32) 

Given the reference domain -D(O) has zero strain and zero stress, (P X y) eq — 0, the change in the free energy will be 
0(j 2 ) and may be ignored in Eq. (J29j) . We now use Eqs. ([20]) . ([30]) & ([29)) to obtain the response to the impulse, 

(P*v(t)) = -7i<<?oo) ,e 9 +^7i<^) , e9 

-p Vll (P xy (-t)P xy (0)) Qeq . (33) 

At time t = + , Eq. ([33"]l coincides with the stress that one would calculate for a sudden impulse (ie the response is 
given by the infinite frequency shear modulus alone). At very long times, t — > oo (where the autocorrelation function 
fully decays), Eq. ([33]) reduces to that given by the zero frequency shear modulus Eq. ([20]) . At first sight, Eq. ([29]) 
looks paradoxical. The first term on the right hand side is the correct long time answer for any changes in the shear 
rate that are completed before time t. Since that first term has no memory it looks as though it gives the infinite 
frequency response. However, as this example has just proved, this is not the case. 

If we apply Eq. (f33[) to a fluid, which must have (goo) eq — flV (P xy ) , then Eq. ([33]) will be compatible with Eq. 

(Ell). 

3. Initially Strained crystal 

For a crystal which is initially strained by an amount 70 = 7(0~), but that is none the less in equilibrium, Eq. lj30j) 
becomes 

(^v> 70+7l ,e 8 = <^t/> 70 , e , " G (7o)7i + 0(7i 2 ), (34) 
and for the change in free energy we have 

—g^~ = - V < P ^)7o,e 9 ■ ( 35 ) 
We now use Eq. ([29]) and consider the response to an impulsive change in the strain rate, j(t) = S(t)-fi. We see that 
(P xy (t)) = (P xy ) ioeq ~G ( 10 h 1 ~/3V 11 (AP xv (T(-t))AP xv (T)) 7oeq 

= (PxvKm ~ ^ (9°°K,e q + Wli ( AP *y) yo , eq - Wli (AP xy (T(-t))AP xy (T)) 7Q eq , (36) 
where AP xy (T) = P xy (T) — (P X y)~ eq - Again we see that this equation gives the initial response to the impulse, 

(P xy (Q + )) = {P xy ) 10 . eq ~ 7i (9oo) 10<eq , (37) 
at time + and as t — * 00 it decays to the equilibrium value, 

(P*y) lQ+ll , eq = (P*y) 70 , eq ~ 7i <5oo) 70 , eg + /3F 7l <Ai£,> • ( 38 ) 



= (^„> 7o , e ,-Go(7o)7i- (39) 

For the case where 70 = we have (P X y) l0 eq = and AP xy {T) = P xy (T) and Eq. ([36]) reduces to Eq. ([33]). We 
reiterate that the values for the various elastic moduli will of course depend on the alignment of the crystal relative 
to the strain rate tensor. For a crystal the elastic modulus is in fact a 4 th rank polar tensor. 



D. Oscillatory Planar Shear 



We now consider the case of oscillatory planar shear applied to a crystal, by using Eq. I|29p to calculate the response 
in the stress B(T) — P xy to an applied strain of the form, 

7 (t) = 70 sinM) = -to {i l0 e luJt } . (40) 

The response to the oscillatory strain will become sinusoidal after the decay of initial transients and is often expressed 
in terms of the storage (or real) Gr and loss (or imaginary) Gj shear moduli 



lim (P, 

t— s-OG 



> xy (t)}=to{i'r G(Lj)e iu ' t } (41) 



where G(u>) = Gr(uj) + iGi{u>), This quantity is related to the complex frequency dependent shear viscosity by the 
equation, 

G(w) = iujfj{uj), (42) 
where fj(u>) — rj^(uj) — irji(w) and thus Go = G(uj — 0) and G^ = lim G(lu). 

uj — >oc 

The applied field is F e = 7(f) = luj cos(wf), the change in free energy for a crystal is given by dA/d-f — — (P xy ) eq = 
and the flux is J(T) = P xy (T). Using Eq. $29$) we have 

(P xy (t)) = {P xy ) l{t)m - /3^ 70 f ds (P ay (-«)Px If (0)) 7(0))efl cos( w (f - a)) (43) 

o 

which, using a trigonometric identity and Eq. (|30|) . gives 

(Px y (t)) = -G 7osm(wt) (44) 

~f3Vuj^f sin(wt) / ds (P xy (-s)P xy (0)) /y{a)eq sm(ujs) 

•J 

-PVuj-fo cos(wf) / ds (P xy (-s)P xy (0)) i{0) e(} cos(ujs). 

Combining this with Eq. (|4l|) we obtain, 

G r (uj) = Go+f3Vuj / ds sm(ujs)C(s) 
Jo 

/>oo 

G/(w) = (3Vuj / ds cos(ws)G(s), (45) 
Jo 

where the correlation function is given by 

C(s) = (P xy (-s)P xy (Q)) eqfi . (46) 
Using Eqs.ldU) and (@5|) we see that, T]r{<jS) — Gi(uj)/lj and Gr(uj) = Go + u)rji((j). 

rjn{uj) = f3V / ds cos(ws)G(s) 



777(0;) = -Gq/w + PV I ds sin(ws)G(s). (47) 
Jo 

If we now ask what memory function 77(f), 

(P xy (t)) = - f ds V (t-s)i(s), (48) 



generates this spectrum for the frequency dependent shear viscosity, we see that it is, 

77(f) = G +0VC(t), f>0; (49) 

77(f) = -G , f < 0. 

A double sided Fourier transform of this function gives Eq. ()47|) . The memory function for the zero frequency elastic 
response must be odd in time because a constant strain is an even function of time while a constant strain rate is of 
course odd. For linear elasticity and linear viscosity to both apply to the system, both the strain and the strain rate 
must be small at all frequencies including near zero. 



III. SIMULATION, RESULTS AND DISCUSSION 



A. Simulation Details 

To test the theory we used both equilibrium time correlation data, obtained using Eqs. ([T]) with F e = and 
nonequilibrium molecular dynamics (NEMD) data, obtained from Eqs. (J4j) from oscillatory strain simulations. The 
simulations used N — 108 particles with periodic boundary conditions to model a perfect crystal at a finite tempera- 
ture. Because we model a single crystal with no defects there will be no long range stresses and no large system size 
effects. We know from the pioneering studies of a realistic potential for argon by Barker et. al. Barker et al.— that 
for perfect crystals, 108 atoms is sufficient for quantitative agreement with experiment. 

An equilibrium face centre cubic (FCC) crystal was formed as a cubic array of periodic molecular dynamics unit 
cells. The FCC crystal is commensurate with the simulation cell and the edge of the cube is in the (1,0,0) direction 
of the crystal. A pairwise additive WCA potential, 

(50) 

was used fo r the interaction between the particles. The energy unit is e, the length unit is a, and the time unit 
is yj ma 2 /e where m is the mass. The system has a number density of p = Na 3 /V = 1.15, which results in an 
equilibrium crystal at the two temperatures studied, T = 0.5 and T = 2.5 (e/ks)- 

For the nonequilibrium simulations eight different frequencies where simulated. For those with a frequency of 
to = 2.513 or higher, the maximum amplitude of the strain (as defined in Eq. (HQ} ) was 70 = 0.025, for the lower 
frequencies 70 = 0.08 was used. As the frequency is lowered, with fixed maximum strain amplitude, the signal to 
noise ratio for the loss component of the shear modulus deteriorates because the strain rate goes towards zero. The 
larger amplitude used at low frequencies helped alleviate this problem a little. 

The equations of motion were integrated using a fourth order Runge-Kutta integrator. A time step of dt = 0.002 was 
used for all simulations except for the NEMD simulations at the highest frequency where a time step of dt — 0.001 
was used. All the simulations, both NEMD and equilibrium, were repeated 1000 times to reduce the statistical 
uncertainties. At the highest frequency the NEMD simulations were given at least 6.3 time units to relax to the 
periodic state and for the lowest frequency this was extended to 2000 time units. For all but the lowest three 
frequencies the NEMD simulations were run for 10 periods to obtain the response. At the lowest frequency the 
NEMD simulations were run for only 2 periods. The longest duration used for producing data from the equilibrium 
simulations was 800 time units. 



B. Simulation Results and Discussion 

1. Equilibrium data 

To calculate the response of the system using Eqs. f45|) or f47j) we need to first determine the equilibrium correlation 
function, Eq. (|4"6*)) . and the zero frequency modulus, Eqs. lfT9l) & l(20|) . using data from equilibrium simulations. The 
equilibrium correlation functions, normalised by the temperature, can be seen in Fig. Q] for the temperatures of 
T = 0.5 and T — 2.5. At zero delay time, t — 0, the height of the the correlation function for the higher temperature 
is approximately a factor of 4 larger than that for the lower temperature. The viscosity is given by the area under 
the curve, lim^^oVR^) = Vr{0 + ) — 0V J °° dsC(s), and as it turns out the values obtained for the two different 
temperatures are very similar. For the temperature of T — 2.5 we have fjn(0 + ) — 283 and for the the temperature 
of T = 0.5 we have 7?ij(0 + ) = 222. Although the area under the curve for the higher temperature appears much 
larger in Fig. [IJ by noting the logarithmic time axis it can be seen that the difference decays very rapidly. This is 
the reason the values are relatively close to each other. The limiting viscosity for the crystal increases, rather weakly, 
with temperature. This is in contrast to a liquid where the viscosity decreases with increasing temperature, often 
strongly, but is similar to a dilute gas where the viscosity also increases with temperature. 

One cannot measure the viscosity of a crystal at zero frequency by subjecting it to an unbounded strain. At zero 
frequency, for large enough strain, a crystal will exhibit a nonlinear response, undergo plastic deformation, cleavage 
or fracture. The viscosity we calculate is the limiting zero frequency shear viscosity, lim(w — ► 0). If we subject the 
crystal to an oscillating strain, the viscosity characterises the dissipation of energy in the low frequency limit. 
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Figure 1: The stress correlation function for the FCC crystal at temperatures of T = 0.5 and T = 2.5. Plotted is the function 
/3C(t) = P {P xy (t)Pxy{0)) . Note that the time axis is logarithmic and that the curve with the higher initial value (close to 
6 x 10 3 ) corresponds to the temperature of T = 2.5. 

An alternative way to measure the limiting zero frequency shear viscosity is to subject a crystal to a fixed but very 
small strain rate for a limited period of time t\ which is inversely proportional to the strain rate, t\ = 7 TO /7- If we se t 
the maximum strain, 7 m to be 7 m <~ 0.1 - according to Lindemann's criterion, then we will not cleave or otherwise 
damage the crystal and we will remain in the linear response regime for both the shear rate which is always very small 
and for elastic deformation. As the strain rate decreases towards zero the amount of time we may strain the crystal 
diverges to infinity. So even though the maximum strain the crystal is limited (7 < 7 m ) as the strain rate is reduced 
we have ever more time available for the shearing system to relax to a nonequilibrium. The shear viscosity of this 
limiting steady state is what we call the limiting shear viscosity of a crystal. Of course as the strain rate is lowered the 
deteriorating signal to noise ratio demands that the size of the system, or the number of times the shearing protocol 
is repeated must be increased. This may not be a practical way to measure the limiting viscosity of a solid. 

2. Frequency Dependent Modulus and Viscosity at T = 0.5 

We examine the case for the FCC crystal at a temperature of T = 0.5 in more detail. The correlation function 
C(t) = (P xy (t)P xy (0)) was calculated from equilibrium simulations as discussed above, and Fourier transformed, 
according to Eqs. (|45|i. using the first order Filon's quadrature given in the Appendix, to obtain the storage and 
loss moduli. The value for Go was obtained from the equilibrium simulations using Eqs. (fl9|) & (|20|) . We then 
performed nonequilibrium simulations at various frequencies, and applied a least squares fit of a sinusoidal function 
to the response allowing us to obtain estimates of G at a small number of distinct frequencies. The results of this 
are shown in Fig. [2j It can be seen that the agreement between the two data sets is very good. This confirms the 
correctness of our theoretical expressions for the frequency dependent elastic moduli. In general the low frequency 
Gj data is difficult (or computationally expensive) to obtain reliably due to the small amplitude of the strain rate. 
As mentioned above the strain amplitude is fixed, and therefore the strain rate becomes very small at low frequencies 
with 7 ~ u>. It can be seen that the low frequency data, from the transformed C(t), decays as lim^^o G[(lo) oc oj 
(this results in a gradient of unity for In [Gj] vs \n[v] at low frequencies, see Fig. [2|). This is obvious from Eq. (pSJl. 
upon taking the small angle approximation cos(ws) = 1 + 0(ui 2 ), we see that lim w ^ Gj{oS) = w?7_r(0 + ). 

In Fig. [3] the same data is shown for the complex frequency dependence of the shear viscosity. The real part of 





Figure 2: The storage Gr(ui) and loss Gi(lu) moduli for the temperature of T = 0.5. The storage modulus Gr(lu) is the 
symbols (+) while the loss modulus is the symbols (x). The solid curves are obtained from linear response theory using 
data obtained from the equilibrium simulations while the symbols were obtained directly from the nonequilibrium molecular 
dynamics simulations. 




Figure 3: Frequency dependent viscosity for the crystal at T = 0.5. The symbols (+) represent the real part of the viscosity 
t)r{ui) and the symbols (x) represent the imaginary part rji(uj). The solid lines are from the linear response theory and the 
symbols are from NEMD simulations. The symbols for the real part at the lowest two frequencies are not very accurate. 



the viscosity, fju(u>), converges to the finite value at zero frequency which is given by the area under the correlation 
function shown in Fig. [TJ What is distinctly different in this graph, relative to equivalent data for a typical fluid, is the 
behaviour of the imaginary part of the viscosity fji(u)). We see that this quantity diverges, apparently to infinity, as the 
frequency approaches zero. This is due to a solid having a nonzero value for the zero frequency shear modulus G(0). 
In a fluid the zero frequency value of the imaginary part of the shear viscosity is zero. The distinctive characteristic 
of the frequency dependent viscosity of a solid phase, relative to a fluid phase, is the contrasting behaviour of the 
imaginary part of the viscosity as the frequency approaches zero. For solids the imaginary part diverges to infinity 
while in fluids it decays to zero. 

IV. CONCLUSIONS 

We have derived a set of theoretical expressions for the linear viscoelastic properties of crystals. Computer sim- 
ulations have been carried out which compare the results of direct nonequilibrium molecular dynamics calculations 
for these properties, with the linear response theory expressions calculated using equilibrium simulation data. The 
agreement between these two sets of results confirms the correctness of our expressions for the linear response. 

A glass is often defined as a supercooled liquid with a shear viscosity that is greater than 10 13 poise Debenedetti 
and Stillinger— . In our units, assuming that our potential gives a very approximate model for argon, this would 
correspond to a viscosity of approximately 10 16 . However we should point out that when this statement is made it 
is also assumed that the shear modulus of the supercooled liquid/glass is zero. Once we enter the glass phase the 
shear modulus is actually nonzero and then according to our theory the Green-Kubo expression for the shear viscosity 
changes reducing the magnitude of the shear viscosity somewhat. 

If the systems we studied here are taken to represent argon we have shown via equilibrium and nonequilibrium 
molecular dynamics calculations that the limiting zero frequency shear viscosity rjn(0 + ) of an argon crystal is only 
two orders of magnitude greater than liquid argon at its triple point. [In the units used in this paper the shear viscosity 
of triple point liquid argon is approximately 3.5.] So in contrast to a glass, the crystalline state has no anomalously 
high shear viscosity. Although we have only calculated this viscosity for a single relative alignment between the crystal 
axes and the strain rate tensor, we do not expect that varying this alignment would lead to an increase in viscosity of 
many orders of magnitude. We know of no other work, experimental or theoretical, that has calculated the limiting 
shear viscosity of a crystal. As far as we are aware, all such work refers to the real and imaginary parts of the shear 
modulus. 

The crystal we have studied is a soft inert gas crystal. It is interesting to speculate about the comparative values of 
the limiting shear viscosities of ionic or covalent crystals. Somewhat counter intuitively these types of crystals could 
have limiting viscosity values that at room temperature, may be lower than that of soft inert gas crystals. What is 
important in increasing the viscosity of crystals is anharmonicity in the nearest neighbour forces. Strong, high Q, 
harmonic crystals can be expected to have low shear viscosities. 

Our work also points out the distinctiveness of glassy systems. Glassy systems by definition exhibit an anomalously 
high viscosity. This is quite different from the behaviour of crystalline solids and of the liquids from which a glass 
may be formed. 

We also see an interesting set of contrasting qualitative behaviour for the temperature dependence of the limiting 
shear viscosity. Gases exhibit a positive temperature coefficient, liquids have a negative temperature coefficient while 
crystals again exhibit a positive temperature coefficient for the limiting shear viscosity. 

The zero frequency shear modulus is profoundly different between a crystal and a fluid. For fluids the shear modulus 
is precisely zero whereas in any solid, including crystals, the shear modulus is nonzero. As Max Born stated in 1939, 
reference [9j, "..there can be no ambiguity in the definition of, or the criterion for, melting. The difference between a 
solid and a liquid is that the solid has elastic resistance to shearing stress while a liquid does not." Our work strongly 
supports this assertion but points out that the contrast in behaviour between fluids and crystals shows its greatest 
effect when one compares the limiting zero frequency component of the imaginary part of the shear viscosity. In fluids 
that value is zero while in crystals the limiting value is infinity! 

Appendix: Transforming the autocorrelation function 

A trapezoidal version of Filon's quadrature was used to transform the correlation functions, 



dsf(s) cos(ws) = — (/(fi)sin(wii) - /(t ) sin(wt )) 
1 /(ti) - /(*o) 



LU 2 tl — tQ 



1 



( cos(orfi) — cos(wto)); (51) 



ds/(s) sin(ws) = -(/(to) cos(urf ) - cos(wti)) 



— to 

Using these quadrature it is straightforward to prove, given /(oo) = 0, that 



1 f{h) /W (gja^) - s in(^o)). (52) 



lim a; / dsf(s) cos(ws) = (53) 



lim w / <2s/(s) sin(ws) = /(0) (54) 
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